#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"
extern double density_above,density_below,density_middle;
extern double compute_heaviside_general(),smooth_front();

extern double b00[imax][jmax],b01[imax][jmax],b02[imax][jmax],
         b10[imax][jmax],b11[imax][jmax],b12[imax][jmax],
         b20[imax][jmax],b21[imax][jmax],b22[imax][jmax];
extern int km,kmx,kmy;

int pwr[11],mxk[11],myk[11],ixlo[11],ixhi[11],
    iylo[11],iyhi[11],sumk[11];
double hk[11],wssor[11];
static logical interior_flag;

#define for_all_xk for (ix=ixlo[k]-1;ix<=ixhi[k]+1;ix++)
#define for_all_yk for (iy=iylo[k]-1;iy<=iyhi[k]+1;iy++)
#define for_all_pointsk for_all_yk for_all_xk

#define for_all_xk_int for (ix=ixlo[k];ix<=ixhi[k];ix++)
#define for_all_yk_int for (iy=iylo[k];iy<=iyhi[k];iy++)
#define for_all_pointsk_int for_all_yk_int for_all_xk_int


#define pointk_x(k,x1) (( (x1)>ixlo[k]-1 ) && ( (x1)<ixhi[k]+1 ))
#define pointk_y(k,y1) (( (y1)>iylo[k]-1 ) && ( (y1)<iyhi[k]+1 ))
#define pointk_xy(k,x1,y1) ((pointk_x(k,x1)) && (pointk_y(k,y1)))

#define pointkx_ok (pointk_x(k,ix))
#define pointky_ok (pointk_y(k,iy))
#define pointk_ok  (pointk_xy(k,ix,iy))

double get_hxk(ix,k)
int ix,k;
{
  if ((ix<ixlo[k])&&(!use_symmetric)&&(!use_axisymmetric))
   return farx+hk[k];
  else if (ix>=ixhi[k])
   return farx+hk[k];
  else
   return hk[k];
}

double get_hyk(iy,k)
int iy,k;
{
  if (iy<iylo[k])
   return fary+hk[k];
  else if (iy>=iyhi[k])
   return fary+hk[k];
  else
   return hk[k];
}

double dvalk(ix,iy,k)
int ix,iy,k;
{
int ixk1,ixk2,iyk1,iyk2;

  if (k==1)
   { ixk1=ixk2=ix;iyk1=iyk2=iy; }
  else
   {
    ixk1=pwr[k-2]+pwr[k-1]*(ix-1);
    ixk2=ixk1+1;
    iyk1=pwr[k-2]+pwr[k-1]*(iy-1);
    iyk2=iyk1+1;
   }

#define dist_avg(gdtype) (gdtype(ixk1,iyk1)+gdtype(ixk1,iyk2)+ \
  gdtype(ixk2,iyk1)+gdtype(ixk2,iyk2))/4  

  return smooth_front(dist_avg(gd1),dist_avg(gd2),dist_avg(gd3),
    density_above,density_below,density_middle,density_spread*pwr[k-1]);
  /* is pwr[k-1] the trick? */
}

double rvalk(ix,k)
int ix,k;
{
int ixk1,ixk2;

  if (k==1)
   { ixk1=ixk2=ix; }
  else
   {
    ixk1=pwr[k-2]+pwr[k-1]*(ix-1);
    ixk2=ixk1+1;
   }
  return ( (gr(ixk1)+gr(ixk2))/2.0 );
}

void construct_a_save(k)
int k;
{
int ix,iy,start_ix,ixcvt;
double d1s,d2s,d3s,d4s,hx0,hx1,hy0,hy1,hx2,hy2,
  hxx0,hxx1,hyy0,hyy1;

  start_ix=sumk[k];
  for_all_pointsk
    {
     ixcvt=ix+start_ix;
#define get_d_val(incx,incy) (dvalk(ix+incx,iy+incy,k))
#define get_r_val(incx) (rvalk(ix+incx,k))

#define coeff1a(x1,x2) (get_d_val(x1,x2)/(4*get_r_val(x1)) )
#define coeff1(x1,x2) (get_d_val(x1,x2)/(2*sqr(hk[k])*get_r_val(x1)) )

     b00[ixcvt][iy]=b01[ixcvt][iy]=b02[ixcvt][iy]=0.0;
     b10[ixcvt][iy]=b11[ixcvt][iy]=b12[ixcvt][iy]=0.0;
     b20[ixcvt][iy]=b21[ixcvt][iy]=b22[ixcvt][iy]=0.0;

/* assume (use_corner_difference) */
     if ((mx_type==FREE_BDRY)||(my_type==FREE_BDRY))
      {
       hx0=hx1=hy0=hy1=hk[k];

       if (mx_type==FREE_BDRY)
        {
         hx0=get_hxk(ix-1,k);
         hx1=get_hxk(ix,k);
        }
       if (my_type==FREE_BDRY)
        {
         hy0=get_hyk(iy-1,k);
         hy1=get_hyk(iy,k);
        }

       hx2=hy2=hk[k];
       hxx0=hx2*hx0;hyy0=hy2*hy0;
       hxx1=hx2*hx1;hyy1=hy2*hy1;

       b00[ixcvt][iy]=d1s=-coeff1a(0,0)*(1/hxx0+1/hyy0);
       b20[ixcvt][iy]=d2s=-coeff1a(1,0)*(1/hxx1+1/hyy0);
       b02[ixcvt][iy]=d3s=-coeff1a(0,1)*(1/hxx0+1/hyy1);
       b22[ixcvt][iy]=d4s=-coeff1a(1,1)*(1/hxx1+1/hyy1);

       b01[ixcvt][iy]=-coeff1a(0,1)*(1/hxx0-1/hyy1)-
                    coeff1a(0,0)*(1/hxx0-1/hyy0);
       b21[ixcvt][iy]=-coeff1a(1,0)*(1/hxx1-1/hyy0)-
                    coeff1a(1,1)*(1/hxx1-1/hyy1);
       b10[ixcvt][iy]=-coeff1a(1,0)*(-1/hxx1+1/hyy0)-
                    coeff1a(0,0)*(-1/hxx0+1/hyy0);
       b12[ixcvt][iy]=-coeff1a(0,1)*(-1/hxx0+1/hyy1)-
                    coeff1a(1,1)*(-1/hxx1+1/hyy1);
      }   /* assume free_bdry */
     else
      {   /* fixed boundary */
       b00[ixcvt][iy]=d1s=-coeff1(0,0);
       b20[ixcvt][iy]=d2s=-coeff1(1,0);
       b02[ixcvt][iy]=d3s=-coeff1(0,1);
       b22[ixcvt][iy]=d4s=-coeff1(1,1);
      }
     b11[ixcvt][iy]=-(d1s+d2s+d3s+d4s);
   }  /* looping all the nodes */

if (FALSE)
 {
  printf("checking matrix k,xlo,xhi,ylo,yhi%5d%5d%5d%5d%5d\n",
    k,ixlo[k],ixhi[k],iylo[k],iyhi[k]);
  for_all_pointsk_int
    {
     ixcvt=ix+start_ix;
     if (b12[ixcvt][iy] != b10[ixcvt][iy+1]) printf("b12 ix,iy%5d%5d\n",ix,iy);
     if (b10[ixcvt][iy] != b12[ixcvt][iy-1]) printf("b10 ix,iy%5d%5d\n",ix,iy);
     if (b22[ixcvt][iy] != b00[ixcvt+1][iy+1]) printf("b22 ix,iy%5d%5d\n",ix,iy);
     if (b00[ixcvt][iy] != b22[ixcvt-1][iy-1]) printf("b00 ix,iy%5d%5d\n",ix,iy);
     if (b02[ixcvt][iy] != b20[ixcvt-1][iy+1]) printf("b02 ix,iy%5d%5d\n",ix,iy);
     if (b20[ixcvt][iy] != b02[ixcvt+1][iy-1]) printf("b20 ix,iy%5d%5d\n",ix,iy);
     if (b01[ixcvt][iy] != b21[ixcvt-1][iy]) printf("b01 ix,iy%5d%5d\n",ix,iy);
     if (b21[ixcvt][iy] != b01[ixcvt+1][iy]) printf("b21 ix,iy%5d%5d\n",ix,iy);

     if (k==1)
     {
      if (b12[ixcvt][iy] != a10[ixcvt][iy+1]) 
       printf("a12 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b12[ixcvt][iy],a10[ixcvt][iy+1]);
      if (b10[ixcvt][iy] != a12[ixcvt][iy-1]) 
       printf("a10 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b10[ixcvt][iy],a12[ixcvt][iy-1]);
      if (b22[ixcvt][iy] != a00[ixcvt+1][iy+1]) 
       printf("a22 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b22[ixcvt][iy],a00[ixcvt+1][iy+1]);
      if (b00[ixcvt][iy] != a22[ixcvt-1][iy-1]) 
       printf("a00 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b00[ixcvt][iy],a22[ixcvt-1][iy-1]);
      if (b02[ixcvt][iy] != a20[ixcvt-1][iy+1]) 
       printf("a02 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b02[ixcvt][iy],a20[ixcvt-1][iy+1]);
      if (b20[ixcvt][iy] != a02[ixcvt+1][iy-1]) 
       printf("a20 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b20[ixcvt][iy],a02[ixcvt+1][iy-1]);
      if (b01[ixcvt][iy] != a21[ixcvt-1][iy]) 
       printf("a01 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b01[ixcvt][iy],a21[ixcvt-1][iy]);
      if (b21[ixcvt][iy] != a01[ixcvt+1][iy]) 
       printf("a21 ix,iy%5d%5d%10.5f%10.5f\n",
         ix,iy,b21[ixcvt][iy],a01[ixcvt+1][iy]);
     }   /* if k==1 */

    }   /* looping check */
  }  /* if false */

 
}   /* construct_a_save */

void zap_v(dest,k)
double dest[imax][jmax];
int k;
{
int ix,iy,start_ix;

  start_ix=sumk[k];
  for_all_pointsk_int
   {
    dest[ix+start_ix][iy]=0.0;
   }
}

void assign_v(dest,source,k)
double source[imax][jmax],dest[imax][jmax];
int k;
{
int ix,iy,start_ix;

  start_ix=sumk[k];
  for_all_pointsk_int
   {
    dest[ix+start_ix][iy]=source[ix+start_ix][iy];
   }

}

void add_v(k,dest,source)
double source[imax][jmax],dest[imax][jmax];
int k;
{
int ix,iy,start_ix;

  start_ix=sumk[k];
  for_all_pointsk_int
   {
    dest[ix+start_ix][iy]+=source[ix+start_ix][iy];
   }

}

logical at_vertexk(ix,k)
int ix,k;
{
  if (
      ((use_symmetric)||(use_axisymmetric)) &&
      (ix<=ixlo[k])  /* test at 1 not 0 due to corner procedure */
     )
   return TRUE;
  else
   return FALSE;
}   /* at_vertexk */

double gstk_general(stk,ix,iy,k)
int ix,iy,k;
double stk[imax][jmax];
{
double m,b;
int ixcvt,start_ix;

  if (pointk_ok)
   return stk[ix+sumk[k]][iy];

  start_ix=sumk[k];
  ixcvt=ix+start_ix;
  if (pointkx_ok)
   {
    if ((my_type==REFLECT_BDRY) || (my_type==SYMMETRIC) ||
        (my_type==FREE_BDRY))
     return 0.0;
    else if (my_type==PERIODIC)
     return ( (iy<iylo[k]) ? stk[ixcvt][iyhi[k]] : stk[ixcvt][iylo[k]]);
    else
     return 0.0;
   }  /*  ix in the square */
  else if (pointky_ok)
   {
    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (at_vertexk(ix,k)) || (mx_type==FREE_BDRY))
     return 0.0;
    else if (mx_type==PERIODIC)
     return ( (ix<ixlo[k]) ? 
       stk[ixhi[k]+start_ix][iy] : stk[ixlo[k]+start_ix][iy]);
    else
     return 0.0;
   }   /* iy in the square */
  else
   {
    if ((mx_type==PERIODIC)&&(my_type==PERIODIC))
     {
      ix=( (ix<ixlo[k]+1) ? ixhi[k] : ixlo[k]); 
      /* ix from ixmin-1 to ixmin above */
      iy=( (iy<iylo[k]+1) ? iyhi[k] : iylo[k]); 
      /* iy from iymin-1 to iymin above */

      return stk[ix+start_ix][iy];
     }
    else
     return 0.0;
   }   /* at a corner */
}    /* gstk_general */

void resid_v(k,uval,fval,rval,error_inf)
int k;
double fval[imax][jmax],uval[imax][jmax],rval[imax][jmax],
  *error_inf;
{
int ix,iy,ixcvt,start_ix;

  start_ix=sumk[k];
  *error_inf=0.0;

  for_all_pointsk
   if (! pointk_ok)
    uval[ix+start_ix][iy]=gstk_general(uval,ix,iy,k);

  /* yloop must be on the outside ! */
  if (k==1)
  {
   for_all_pointsk_int
    {
     rval[ix][iy]=(fval[ix][iy]-
       a12[ix][iy-1]*uval[ix][iy-1]-
       a21[ix-1][iy]*uval[ix-1][iy]-
       a12[ix][iy]*uval[ix][iy+1]-
       a21[ix][iy]*uval[ix+1][iy]-
       a22[ix-1][iy-1]*uval[ix-1][iy-1]-
       a02[ix+1][iy-1]*uval[ix+1][iy-1]-
       a22[ix][iy]*uval[ix+1][iy+1]-
       a02[ix][iy]*uval[ix-1][iy+1]-
       a11[ix][iy]*uval[ix][iy] );
     *error_inf=max(*error_inf,fabs(rval[ix][iy]));
    }  /* looping stencil */
  }  /* k==1 */
  else
  {  /* k>1 */
   for_all_pointsk_int
    {
     ixcvt=ix+start_ix;
     rval[ixcvt][iy]=(fval[ixcvt][iy]-
       b12[ixcvt][iy-1]*uval[ixcvt][iy-1]-
       b21[ixcvt-1][iy]*uval[ixcvt-1][iy]-
       b12[ixcvt][iy]*uval[ixcvt][iy+1]-
       b21[ixcvt][iy]*uval[ixcvt+1][iy]-
       b22[ixcvt-1][iy-1]*uval[ixcvt-1][iy-1]-
       b02[ixcvt+1][iy-1]*uval[ixcvt+1][iy-1]-
       b22[ixcvt][iy]*uval[ixcvt+1][iy+1]-
       b02[ixcvt][iy]*uval[ixcvt-1][iy+1]-
       b11[ixcvt][iy]*uval[ixcvt][iy] );

     *error_inf=max(*error_inf,fabs(rval[ixcvt][iy]));
    }  /* looping stencil */
  }   /* k>1 case */ 
}  /* resid_v */

void gs_iterate(k,fval,uval)
int k;
double fval[imax][jmax],uval[imax][jmax];
{
int ix,iy,ixcvt,start_ix;

  start_ix=sumk[k];
  for_all_pointsk
   if (! pointk_ok)
    uval[ix+start_ix][iy]=gstk_general(uval,ix,iy,k);

  /* yloop must be on the outside ! */
  if (k==1)
   {
    for_all_pointsk_int
     uval[ix][iy]=wssor[k]*(fval[ix][iy]-
      a12[ix][iy-1]*uval[ix][iy-1]-
      a21[ix-1][iy]*uval[ix-1][iy]-
      a12[ix][iy]*uval[ix][iy+1]-
      a21[ix][iy]*uval[ix+1][iy]-
      a22[ix-1][iy-1]*uval[ix-1][iy-1]-
      a02[ix+1][iy-1]*uval[ix+1][iy-1]-
      a22[ix][iy]*uval[ix+1][iy+1]-
      a02[ix][iy]*uval[ix-1][iy+1])/a11[ix][iy]+
       (1-wssor[k])*uval[ix][iy];
    for (iy=iyhi[k];iy>=iylo[k];iy--)
    for (ix=ixhi[k];ix>=ixlo[k];ix--)
     uval[ix][iy]=wssor[k]*(fval[ix][iy]-
      a12[ix][iy-1]*uval[ix][iy-1]-
      a21[ix-1][iy]*uval[ix-1][iy]-
      a12[ix][iy]*uval[ix][iy+1]-
      a21[ix][iy]*uval[ix+1][iy]-
      a22[ix-1][iy-1]*uval[ix-1][iy-1]-
      a02[ix+1][iy-1]*uval[ix+1][iy-1]-
      a22[ix][iy]*uval[ix+1][iy+1]-
      a02[ix][iy]*uval[ix-1][iy+1])/a11[ix][iy]+
       (1-wssor[k])*uval[ix][iy];
   }  /* k==1 */
  else
   {
    for_all_pointsk_int
    {
     ixcvt=ix+start_ix;
     uval[ixcvt][iy]=wssor[k]*(fval[ixcvt][iy]-
      b12[ixcvt][iy-1]*uval[ixcvt][iy-1]-
      b21[ixcvt-1][iy]*uval[ixcvt-1][iy]-
      b12[ixcvt][iy]*uval[ixcvt][iy+1]-
      b21[ixcvt][iy]*uval[ixcvt+1][iy]-
      b22[ixcvt-1][iy-1]*uval[ixcvt-1][iy-1]-
      b02[ixcvt+1][iy-1]*uval[ixcvt+1][iy-1]-
      b22[ixcvt][iy]*uval[ixcvt+1][iy+1]-
      b02[ixcvt][iy]*uval[ixcvt-1][iy+1])/b11[ixcvt][iy]+
       (1-wssor[k])*uval[ixcvt][iy];
    }
    for (iy=iyhi[k];iy>=iylo[k];iy--)
    for (ix=ixhi[k];ix>=ixlo[k];ix--)
    {
     ixcvt=ix+start_ix;
     uval[ixcvt][iy]=wssor[k]*(fval[ixcvt][iy]-
      b12[ixcvt][iy-1]*uval[ixcvt][iy-1]-
      b21[ixcvt-1][iy]*uval[ixcvt-1][iy]-
      b12[ixcvt][iy]*uval[ixcvt][iy+1]-
      b21[ixcvt][iy]*uval[ixcvt+1][iy]-
      b22[ixcvt-1][iy-1]*uval[ixcvt-1][iy-1]-
      b02[ixcvt+1][iy-1]*uval[ixcvt+1][iy-1]-
      b22[ixcvt][iy]*uval[ixcvt+1][iy+1]-
      b02[ixcvt][iy]*uval[ixcvt-1][iy+1])/b11[ixcvt][iy]+
       (1-wssor[k])*uval[ixcvt][iy];
     }
    }  /* k>1 */

}  /* gs_iterate */

void coarse_v(k,uval,uvalcoarse)
int k;
double uval[imax][jmax],uvalcoarse[imax][jmax];
{
int ix,iy,start_ix,start_ix_coarse,ixcvt,iycvt;

  start_ix=sumk[k];
  for_all_pointsk
   if (! pointk_ok)
    uval[ix+start_ix][iy]=gstk_general(uval,ix,iy,k);
  k++;
  start_ix_coarse=sumk[k];
  for_all_pointsk_int
   {
    ixcvt=2*ix+start_ix;
    iycvt=2*iy;
    uvalcoarse[ix+start_ix_coarse][iy]=(uval[ixcvt-1][iycvt-1]+
      uval[ixcvt-1][iycvt+1]+uval[ixcvt+1][iycvt-1]+
      uval[ixcvt+1][iycvt+1]+2.0*(uval[ixcvt+1][iycvt]+
      uval[ixcvt][iycvt-1]+uval[ixcvt-1][iycvt]+
      uval[ixcvt][iycvt+1])+4*uval[ixcvt][iycvt])/16.0;
   }
}   /* coarse_v */

void fine_v(k,uval,uvalfine)
int k;
double uval[imax][jmax],uvalfine[imax][jmax];
{
int ix,iy,ixcvt,iycvt,start_ix,start_ix_fine,ixcvt1;

  start_ix=sumk[k];
  for_all_pointsk
   if (! pointk_ok)
    uval[ix+start_ix][iy]=gstk_general(uval,ix,iy,k);
  start_ix_fine=sumk[k-1];

  for_all_pointsk   /* MUST include exterior coarse pts for interior fine */
    {
     ixcvt=2*ix+start_ix_fine;
     iycvt=2*iy;
     ixcvt1=ix+start_ix;
     uvalfine[ixcvt][iycvt]=uval[ixcvt1][iy];
     uvalfine[ixcvt+1][iycvt]=(uval[ixcvt1][iy]+uval[ixcvt1+1][iy])/2.0;
     uvalfine[ixcvt+1][iycvt+1]=(uval[ixcvt1][iy]+uval[ixcvt1+1][iy]+
      uval[ixcvt1][iy+1]+uval[ixcvt1+1][iy+1])/4.0;
     uvalfine[ixcvt][iycvt+1]=(uval[ixcvt1][iy]+uval[ixcvt1][iy+1])/2.0;
    }
}     /* fine_v */

void initgrids_v()
{
int k;

 km=min(kmx,kmy);

 pwr[0]=1;
 sumk[0]=sumk[1]=sumk[2]=0;
 for (k=1;k<=km;k++)
  {
   pwr[k]=2*pwr[k-1];
   mxk[k]=mx/pwr[k-1];
   myk[k]=my/pwr[k-1];
   hk[k]=h*pwr[k-1];
   ixlo[k]=iylo[k]=1;
   ixhi[k]=mxk[k]-1;
   if (ixmax==mx) ixhi[k]++;
   iyhi[k]=myk[k]-1;
   if (iymax==my) iyhi[k]++;
   if (k>2)
    sumk[k]=sumk[k-1]+(ixhi[k-1]-ixlo[k-1]+3);
   if (k>1)
    construct_a_save(k);
   wssor[k]=2.0/(1.0+2.0*sin(pi*0.5/min(mxk[k],myk[k])));
   wssor[k]=1.0;
  }   /* looping all k */
}    /* initgrids */
